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The transfer function T{k) of dark matter (DM) perturbations during matter domination is 

QQ ■ obtained by solving the linearized collisionless Boltzmann-Vlasov equation. We provide an exact 

f^ ' expression for T{k) for arbitrary distribution functions of decoupled particles and initial conditions, 

f^ , which can be systematically expanded in a Fredholm series. An exhaustive numerical study of 

CNj . thermal relics for different initial conditions reveals that the first two terms in the expansion of 

jp». ' T{k) provide a remarkably accurate and simple approximation valid on all scales of cosmological 

fl\ [ relevance for structure formation in the linear regime. The natural scale of suppression is the 

^^ ' free streaming wavevector at matter-radiation equality, fc/s(feq) = l^'n'PoM /[{V'^) (1 + Zeq)]\ . An 

important ingredient is a non-local kernel determined by the distribution functions of the decoupled 
particles which describes the memory of the initial conditions and gravitational clustering and yields 
a correction to the fluid description. This correction is negligible at large scales k <^ kfs{teq) but 
it becomes important at small scales k > kfs{teq). Distribution functions that favor the small 
momentum region yield longer- range memory kernels and lead to an enhancement of power at small 
scales k > kf^it^q). Fermi-Dirac and Bose-Einstein statistics lead to long-range memory kernels, 
Q ' with longer range for bosons, both resulting in enhancement of T{k) at small scales. For DM thermal 

,y , relics that decoupled while ultrarelativistic we find kfsiteq) — 0.003 {gd/'2)'^ (m/keV) [kpc]~ , 

^ ■ where gd is the number of degrees of freedom at decoupling. For WIMPS we obtain kfs{teq) = 

5.88 (5d/2 )3 (m/100GeV)5 (T^/lOMeV)^ [pc]-\ For k < kj.it^q), T{k) ~ 1 - C[k/kf,{t,q)]^ 

where C = 0(1) and independent of statistics for thermal relics. We provide simple and accurate fits 

04 ' for T{k) in a wide range of small scales k > kfg{teq) for thermal relics and different initial conditions. 

^ ' The numerical and analytic results for arbitrary distribution functions and initial conditions allow 

Cn , an assessment of DM candidates through their impact on structure formation. 
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I. INTRODUCTION AND RESULTS 

The concordance ACDM standard cosmological model successfully explains a wide range of highly precise astro- 
physical and cosmological observations. Main ingredients are an early stage of accelerated expansion (inflation) a 
more recent stage of accelerated expansion driven by dark energy, and the presence of dark matterfDM) composed 
of primordial particles which are cold and coUisionless at the time when the first structures formed[l|,[4|. 

In the cold dark matter (CDM) model, structure formation proceeds in a hierarchical bottom up approach: small 
scales become non-linear and collapse first and their merger and accretion leads to structure on larger scales. This is 
a consequence of the fact that CDM features negligible small velocity dispersions leading to a power spectrum that 
favors small scales. In this hierarchical scenario, dense clumps that survive the merger process form satellite galaxies. 

Large scale ACDM simulations seemingly lead to an overprediction of satellite galaxies [3| by almost an order 
of magnitude over the number of satellites that have been observed in Milky- Way sized galaxies [3, 0, H, B 0]- 
These simulations also yield a distinct prediction: virialized DM halos should feature a density profile that increases 
monotonically towards the center [J, [^, [3, \1Q, [ll| such as the Navarro- Frenk- White (NFW) profile [a] or more general 
central density profiles p{r) ^ r~^ with 1 < /3 ^ 1.5 [a, SIlll- These profiles accurately describe clusters of galaxies 
but indicate a divergent cusp at the center of the halo. 

There is, however, an accumulating body of observational evidence [l^, [iJ, Hj U3, Ss [l3, SM \1m that seem to 
indicate that the density profile in the central regions of dwarf galaxies is smooth leading to the suggestion that the 
central regions feature smooth cores instead of cusps as predicted by CDM. 

Warm dark matter (WDM) particles were invoked [20, [2l[, [23] as possible solutions to the above mentioned dis- 
crepancies both in the overabundance of satellite galaxies and as a mechanism to smooth out the cusped density 
profiles predicted by CDM simulations into the cored profiles seemingly observed in dwarf-spheroidal satellite galaxies 
(dShps). WDM particles feature a range of velocity dispersion in between the CDM and hot dark matter (HDM) 
leading to free streaming scales that smoothes out small scale features and could be consistent with core radii of the 
dSphs. 

Even if these problems of the CDM model find other astrophysical solutions, it remains an important and intrinsically 
interesting question to assess the clustering properties of WDM candidates, since these may emerge from extensions 
beyond the standard model of particle physics. Furthermore, it remains an important aspect to understand possible 
departures of the CDM paradigm at small scales provided by alternative particle physics candidates. 

The gravitational clustering properties of coUisionless DM in the linear regime are described by the power spectrum 
of gravitational perturbations and in particular the transfer function that converts the primordial power spectrum 
into the power spectrum later in the matter dominated era, which in turn, is the input in the large scale numerical 
simulations of galaxy formation. Free streaming [2| of coUisionless DM leads to a suppression of the transfer function 
on length scales smaller than the free streaming scale via Landau damping f23|. 

Perturbations in a coUisionless system of particles with gravitational interactions is fundamentally different from 
fluid perturbations in the presence of gravity. The (perfect) fluid equations correspond to the limit of vanishing mean 
free path. In a gravitating fluid, pressure gradients tend to restore hydrostatic equilibrium with the speed of sound 
in the medium, and short wavelength fluctuations are simple acoustic waves. For large wavelengths, the propagation 
of pressure waves cannot halt gravitational collapse on a dynamical time scale. The dividing line is the Jeans length: 
perturbations with shorter wavelengths oscillate as sound waves, while perturbations with longer wavelengths undergo 
gravitational collapse. 

In a gas of coUisionless particles with gravitational interaction the situation is different since the mean free path is 
much larger than the size of the system (Hubble radius) and the fluid description is not valid. Instead, the coUisionless 
Boltzmann-Vlasov equation for the distribution function must be solved to obtain the dynamics of perturbations 



Just as in the case of plasma physics, the Hnearized Boltzmann-Vlasov equation describes collective excitations 
[201 . In the case of a colhsionless gas with gravitational interactions these collective excitations describe particles free- 
streaming in and out of the gravitational potential wells of which they are the source. The damping of short wavelength 
collective excitations is akin to Landau damping in plasmas [26[: it is a result of dephasing via phase mixing when 
the particles are out of phase with the potential wells that they produce [23| and leads to the colhsionless damping of 
the collective modes |28| |. 

Gilbert [29| studied the linearized Boltzmann-Vlasov equation in a matter dominated cosmology for non-relativistic 
particles described by an (unperturbed) Maxwell-Boltzmann distribution function. This equation was cast as a 
Volterra integral equation whose numerical solution revealed a limiting value of the wavevector above which density 
perturbations are Landau damped and below which they grow via gravitational (Jeans) instability [29]. The results 
were consistent with replacing the speed of sound by the Maxwellian velocity dispersion in the Jeans length [up to a 
normalization factor of 0(1)]. Gilbert's equations were solved numerically to study: (i) the colhsionless damping of 
density fluctuations in an expanding cosmology with massive neutrinos with an approximated Fermi-Dirac distribution 
function [23| ; (ii) the dissipationless clustering of neutrinos with a truncation of the Fermi-Dirac distribution function 
and an analytic fit to the numerical solution of the integral equation [33|- These results were also used to analyze 
the linear regime [21i[ ; (iii) cosmological perturbations and massive light neutrinos [25|, l3l| , and more recently similar 
integral equations were solved approximately for thermal neutrino relics [32l [. 

DM may be composed of several species [33], a possibility that can be accommodated in most extensions of the 
Standard Model, with candidates ranging from weakly interacting massive particles (WIMPs) [3] to axions [2] and 
sterile neutrinos [3J, [Sa, [30] . DM candidates may be produced in the early Universe via different mechanisms and 
some of them are conjectured to decouple and freeze-out with non-thermal distribution functions. Sterile neutrinos 
are an example of this DM candidate [3j, [35|, [36| . It is important to understand the clustering properties of possible 
DM particles that decoupled in or out of local thermal equilibrium. Predictions for the effective free-streaming length 
and transfer function are a necessary ingredient in the study of structure formation. The free-streaming length in 
such mixture was recently investigated in ref. [37] for arbitrary distribution functions of relic particles that decoupled 
in or out local thermodynamic equilibrium (LTE). This study analyzed the Boltzmann-Vlasov equation in Minkowski 
space-time focusing on the marginally unstable collective modes. It revealed important (and intriguing) aspects of 
the free streaming lengths associated with the distribution function of the decoupled particles. 

A program that bridges the microphysics of the production, evolution and decoupling of DM particles with large 
scale structure formation begins by obtaining the distribution function of the DM particles from the dynamics of 
production and decoupling. The DM distribution function determines the abundance [2], primordial phase space 



properties and generalized constraints on the DM candidate [38|, [39| and their free streaming lengths |37| . Results 
in this direction have been reported [39l . [40| . The primordial phase space density of particles T) permits to obtain 
deep insights on D M [39l . The phase space density can only decrease during gravitational dynamics via mergers and 
"violent relaxation" [4l| . Recent photometric and kinematic data on dwarf spheroidal satellite galaxies in the Milky 
Way (dShps) and the observed DM density today yield upper and lower bounds on the mass, primordial phase space 
densities and velocity dispersion of the DM candidates 39* . Combining these constraints with recent results from 
A^-body simulations yield estimates on the mass of the DM particles in the range of a few keV ^39] . 

The DM distribution function after freeze-out is the main ingredient to obtain the transfer function and the power 
spectrum in the linear theory. In prin ciple the DM transfer function may be obtained by modifying the publicly 
available CMB anisotropy codes[42,[43] that treat baryons, photons and dark matter [15], to account for the different 
distribution functions of several species and range of parameters, masses and couplings. 

In practice, this is a computationally daunting problem: fairly complicated non-equilibrium distribution functions 
for a variety of possible DM species (axions, sterile neutrinos, etc) must be input in the codes (with substantial 
modifications of the standard codes). The exploration of the parameter space (masses, couplings and mixing angles) 
and their impact on the transfer function would require an enormous computational effort. 

After matter-radiation equality, the contribution from the baryon-photon fluid to the DM transfer function can only 
be a few percent[l|. DM only couples to the baryon-photon fluid via gravitational interaction and the most prominent 
feature of this coupling on the DM transfer function arc baryon acoustic oscillations (BAO) on a scale corresponding 
to the sound horizon at recombination, corresponding today to about 150Mpc[45, J^, 47]. Our goal is to understand 
the (DM) transfer function at much smaller scales A < Mpc, at which the effect of the acoustic oscillations of the 
baryon-photon fluid (BAO) can be safely neglected. 

Objectives: The main goals of this article are the following: 

(a): to provide an analytic, accurate and simple framework to obtain the DM transfer function during matter 
domination for general initial conditions and arbitrary distribution functions of relic DM particles that decoupled in 
or out of LTE and that are non-relativistic after matter-radiation equality. Neglecting the gravitational coupling to 



the baryon-photon fluid entails that the DM transfer function will eventually require corrections at the few percent 
level. However, we seek to gain understanding of robust features of T{k) at small scales. 

(b): to understand the impact of the statistics of the relic particle on the small scale properties of the transfer 
function. More precisely, we study which features of the distribution function affect more prominently the power 
spectrum at small scales. 

Summary of Results: We study the linearized coUisionless Boltzmann-Vlasov equation in the non-relativistic 
limit for particles that decoupled in or out of LTE with arbitrary (but isotropic) distribution functions and general 
initial conditions. 

It is transformed into an integro-differential equation for density perturbations which features non-local kernels that 
include memory of gravitational clustering and represent corrections to the fluid description. 

The influence of this memory becomes more important at small scales. Fermi-Dirac (FD) and Bose-Einstein (BE) 
statistics result in long range memory, with the longer range for Bose-Einstein statistics. Distribution functions that 
favor small momenta lead to longer range memory and result in an enhancement of the transfer function at small 
scales. 

We obtain an exact expression for T(fc) for arbitrary initial conditions and distribution functions and provide a 
Fredholm expansion that systematically includes the effects of memory of gravitational clustering and higher moments 
of the distribution function, which are corrections to the fluid approximation. 

An exhaustive numerical study of thermal relics that decoupled relativistically, (Fermions (FD) and Bosons, (BE)) or 
non-relativistically [WIMPs, (MB)] for different initial conditions reveal that the first two terms of the Fredholm series, 
given by eqn. (|2.11ip provide an accurate and simple approximation to the exact transfer function. These include 
explicitly the corrections to the fluid description and the influence of statistics on small scales. This approximation 
furnishes a useful tool to extract the broad properties of the transfer functions for arbitrary initial conditions and 
distribution functions. 



The approximate expression for the transfer function given by ea. (|2.11ip features three main ingredients: (i): 
the solutions of Jeans' fluid equations, (ii:) the free streaming solution of the Boltzmann equation in absence of 
gravitational perturbations, (iii:) a non-local kernel that depends on the distribution function of the decoupled 
particles, it defines the second contribution in eq. (|2.11ip . includes memory of gravitational clustering and represents 
a correction beyond the fiuid approximation. 

The scale of suppression of T(fc) is the comoving free streaming wave- vector at matter-radiation equality, kfs(teq), 
where 
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and the angular brackets refer to the average with the distribution function of the decoupled particle, poM is the DM 
matter density today and a{t) the scale factor. 
We find, 



i^fsyte 



5.E 



pc 
0.00284 



V 2 / VlOOGeV/ VlQMeV / 



WIMPs 



[kpc] 



lOMeV 

FD thermal relics 



(1.2) 



BE thermal relics 



^.lOOGeV^ 
2 J ^ keV 

l»"''^(f)'i^^- 

where g^ is the number of relativistic species at decoupling. 

For large scales k < kfs{teq), the contribution of the non-local memory kernel is subleading, whereas for scales 
k > kfs{teq) the memory of gravitational clustering and corrections beyond the fluid approximation become important. 

FD and BE distribution functions lead to long range memory kernels with Bose-Einstein statistics leading to the 
longest range. Longer range memory leads to an enhancement of the transfer function (and power spectrum) at small 
scales as depicted in flg. [TOl For k > kfs{teq), (7 > \/2) the difference in statistics of the decoupled particles becomes 
very important 

At large scales k <^ kfs{teq) the transfer function T{k) in terms of the variable 7 — \f2 k/kfs{teq) has the following 
behaviour for all cases considered 

n^) = i-(:^) +0(7') (1-3) 



where 70 is the same for Fermi-Dirac, Bose-Einstein and MaxweU-Boltzmann distribution functions, and is given by 

eq.dMID. 

The transfer function for the three different cases (MB, FD and BE) can be compared by parametrizing them in 
terms of 7. This comparison for different initial conditions is displayed in fig. 1101 Simple functional fits to T{k) are 
provided on different regions of k for various initial conditions for thermal relics. 

II. GILBERT'S EQUATION 

After decoupling from the plasma, or freeze-out, the distribution function of the decoupled particles obeys the 
coUisionless Boltzmann-equation. In absence of gravitational perturbations the solution of this equation yields a 
distribution function of the form [2, [33] 

/o(P/(i),i)=/oUw^U/o(p), (2.1) 

where Pf{t), a{t), ad are the physical momentum, scale factor and its value at decoupling respectively, and p is the 
comoving momentum. Consistently with isotropy we assume that /o is a function of p = \p\. 

During the matter dominated era, for modes that are deep inside the Hubble radius and DM particles that are 
non-relativistic the evolution of DM density perturbations and gravitational perturbations is studied with the non- 
relativistic Boltzmann-Vlasov equation^]- The non-relativistic case applies to all comoving scales (physical scales 
today) A < 170 Mpc which were inside the Hubble radius at matter-radiation equality and DM particles with masses 
m > few eV that decoupled before matter-radiation equality. 

To linear order in perturbations the distribution function of the decoupled particle and the Newtonian gravitational 
potential are P, [H, ^^ 

fip:x;t) = fo{p) + Fi{p;x;t) (2.2) 

f{x,t) = ipo{x,t) + ipi{x,t) , (2.3) 

where /o(p) is the unperturbed distribution function of the decoupled particle, (poix, t) is the background gravitational 
potential that determines the homogeneous and isotropic Friedmann-Robertson- Walker metric and p, x are comoving 
variables. 

The reader is referred to refs. [l|, [23|, [2^, [29|, [30, [Ml for details on the linearization of the non-relativistic coUisionless 
Boltzmann-Vlasov equation. 

In conformal time r and in terms of comoving variables p, x the linearized Boltzmann-Vlasov equation is |23l . |29| . [30| 

1 ■ ^ V^Fi - m Vj^i • Vjj/o = , (2.4) 



along with Poisson's equation, 

AnGm r (Pp 
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It is convenient [29|, [30| to introduce a new time variable s related to the conformal time r by 



^1^1 = -^^^ i 7K5r,Fi{x.r). (2.5) 



ds^ — , (2.6) 

a 

and to take spatial Fourier transforms of the gravitational potential (pi {x, t) and the perturbation Fi [x, r) to obtain 

^^'^l^P'^'^ +±lF,{k,p;s) - ik-^pUp) a\s) v,{k,s) = (2.7) 
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Integrating from the initial time s = s^ to s, the Boltzmann-Vlasov equation (|2.7p yields 
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Fi(fc,p;s) = Pi(fc,p;s,)e-*^(^-'^')+zmfc.Vjr/o(p) / ds' e'^^ ^' ^'^^ a\s') ^^{k,s') . (2.9) 



This equation can be turned into an integral equation for the gravitational potential by multiplying both sides by 
—AttG m/[k'^ a(s)]; integrating in p, and using the relation (|2.8p . We obtain 
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(2.10) 
The right hand side of this equation, the inhomogeneity, is determined by the first term in the solution eq. p.Qp 
and describes the free streaming solution of the Boltzmann-equation in absence of gravitational perturbations. This 
can be seen directly by performing the inverse Fourier transform in k to spatial coordinates, which yields for the 
inhomogeneity the following form: 

T[x - —{s - s.^)] where T[x] = -^ j e"^ "^ Fi{k,p]Si) -— . (2.11) 
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For a species that has become non-relativistic: p/m = v where v is the comoving velocity. Using the relation ea. (|2.6p 
and dr = dt/a(t) we see that 

l-is-s.)^l[^^^dt'^lMv;t) (2.12) 

where v(t) = v/a{t) is the physical velocity of the non-relativistic particle and lfs{v] t) is the comoving free streaming 
distance that the particle with comoving velocity v has traveled during (comoving) time from ti until t [2\. The 
physical distance is obtained by multiplying the above expression by a{t). Therefore the inhomogeneity eq. (|2.1ip in 
eq. (|2.10p is of the form 

v^ dv 

Y^ y^[x-lfs(v;t}\ . 

The variable s, related to conformal time as in ea. (|2.6p . obeys the differential equation 

'^' ' (2.13) 



da a^ H[a) 
where a is the scale factor and H the Hubble parameter. For matter dominated cosmology 

H{a) = ^ , (2.14) 

where 

Hqm = — o — Pom = Hq Qm , (2-15) 

with Hq = 100 h kmsec""'^ Mpc~ is the Hubble parameter today and poM is the matter density today with fim = 0.233 
for DM. 

Normalizing the scale factor to unity today, namely a(0) = 1 and taking Sj = as initial value for the time variable 
s, ea. (|2.13p is integrated to yield 
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The initial value of the scale factor a.i corresponds to the value at matter-radiation equality a.i = a^q = 1/(1 + z^q) 

time variable u as 

2u 



with Zeq — 3050. It is convenient to introduce a dimensionless time variable u as 
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SO that 



u= 1- (^V = 1- J^-i^ , 0< M< l-a|„~ 0.982 (2.18) 
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and the scale factor in terms of the variable u is 
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A. The free streaming length. 



We note that during matter domination, the conioving free streaming distance is given by 
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We see that the second term in lfs{t) scales as 1/ \/a{t). This observation will be important below to establish the 
redshift dependence of the comoving free streaming wave- vector. 

Motivated by the expression for the free-streaming distance ea. (|2.12l) and (|2.20p and by the usual Jeans' wavevector 
for fluids J24| , we introduce the physical free streaming wavevector, 
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The comoving free streaming wavevector is defined as 
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This definition of the free streaming wavevector is analogous to that of the Jeans' wavevector in a fluid replacing the 



(2.23) 



speed of sound by y (V"^) (see the discussion below). 

It is clear from this expression that the comoving free streaming wave vector scales as 
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where, using eq. (|2.22p . 
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The comoving free streaming wavelength is given by: 
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There is a simple interpretation of kfs(ti): consider the comoving free streaming distance that a particle with comoving 
velocity v = y {V'^) travels from the initial time ti = teq until today. Setting s = in eq. (|2.20p and neglecting the 
second term since 1/ yJaiQ) = 1 <C Xj Ja^ ~ 55.2, we find 
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We note that the free streaming wavevector at ie? is related to the free streaming distance that the particle has 
traveled between matter-radiation equality and today. 
The matter density today is ^ 
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^ We consider a species with a single internal degree of freedom, a different value g can easily be included in the final expressions. 



It also proves convenient to introduce the density perturbation 



related to the gravitational potential as 



Mk,s)=m I -0^Fiik,p;s), (2.29) 



^,(k;s)ais)^-^^A{k;s), (2.30) 



and the gravitational potential normalized at the initial time 

mu)^^. (2.31) 

where s{u) is given by ea. (|2.17p . 

The unperturbed distribution function /o (p) is a dimensionless function, and as such it can be written as a function of 
the ratio of the comoving momentum p and the value of the decoupling temperature today T^.q and other dimensionless 
quantities, such as the ratio of the mass of the particle to the decoupling temperature, dimensionless coupling constants 
from the microphysics of decoupling etc., [2, l39|, namely 

Mp) = fo{y;xi,x2---) , y = 7fr- , ^1 = 77^ — ' (2-32) 

where T^ is the decoupling temperature and Tdo ~ Td ad is its value today. Using entropy conservation it follows that 
i 
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where Tcmb — 2.348 x 10 ^ eV is the CMB temperature today and gd is the effective number of relativistic degrees of 
freedom at decoupling. For these distributions 
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Integrating by parts the momentum integral in the left hand side of eq. (|2.10p . integrating over the angle k ■ p = cos 9 
and dividing both sides of the integral equation (|2.10p by the initial value 
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the integral equation (|2.10p becomes 



$(fc, u) - - (1 - uf f n[a (u - u')] ^^^^ du' = (1 - uf I[a u] . (2.35) 



where 



Using eq. (|2.33p and the value 



for non-baryonic dark matter, a becomes 
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Um h^ = 0.105 (2.37) 
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V^d/ TO 



The cosmologically relevant range where this approach appHes goes from scales well inside the Hubble radius (where 
the non-relativistic approximation is valid) till the smallest scales where the linearized approximations are valid 



(lOOOMpc)"^ < fc < (0.01 Mpc)"^ . 

Therefore, the range of the dimensionless variable a results as 
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As it is customary, we only consider distributions spherically symmetric on k. The non-local kernel in eqn. (j2.35p is 
given by 
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where we introduced the normalized distribution function 
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where we have assumed that Fi does not depend on the direction of p. 
The inhomogeneity /[au] obeys the initial conditions 



I[au — 0] — l , -—I[au] 
du 

The density perturbation normalized at the initial time namely 

A(fc;s) 



= 



M = 



5{k ,u) 



A(fc;0) 



is related to $(fc; u) [see ea. (P^ and eg. ipTTg)) ] by 



<i>(fc; u) 



a{u) 



S{k,u) = (l-u)^ d{k,u) 



(2.41) 



(2.42) 



(2.43) 



(2.44) 



(2.45) 



(2.46) 



(2.47) 



Then, from eq. ()2.35|) 6{k, u) obeys the Volterra equation of the second kind 

<5(fc, u)-- I nfa [u - u')] f''^' "2 du' = I\a u] . 
a Jo [1 - ""'J 

which is Gilbert's equation f29s|. 

The initial conditions of the inhomogeneity ea. (|2.44p lead to the following initial conditions on the gravitational 
potential and density perturbations 



S{k,u = 0)^1 , — 5{k,u) 
du 

<P(k,u = 0) = 1 , — $(fc,u) 
du 



u=0 



= -2 



n=0 



(2.48) 
(2.49) 
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B. Exact results from Gilbert's equation 

The long wavelength limit a -^ of eq. (|2.47p affords an exact solution. In this limit, 

lim -U\a(u-u')]=u-u' , /[O] = 1 , (2.50) 

Q— >o a 

and {u—u') Q{u—u') is the Green's function of the differential operator d? /d?u. Therefore taking the second derivative 
with respect to u of the Volterra equation (|2.47p , we obtain 

5{0.u)-—-^5{Q,u)^Q, (2.51) 

(1 — uy 

the solution satisfying the initial conditions cq. (|2.48p is 

.(0,.) = ^^ + ^(l-.)3, (2.52) 



which in terms of the scale factor eq. (|2.19p is given by 



5(0,.). I ^ + ^ ^ ^. (2.53) 



O d^q O 



a{u) 



This is the usual long-wavelength solution for density perturbations in a matter dominated cosmology. From the 
eq. (|2.46p the long- wavelength limit of the gravitational potential is 

*i>(0,u) = ^ + ^(l-u)5. (2.54) 

We now show that 5{k, u) behaves as 1/(1 — w)^ for all values of a when u ^ 1. 
Let us assume a general power law behavior for 6lk, u) as u ^ 1, 

S{k, u) "S^ Aia) (1 - u)^'^ [1 + 0(1 - u)] . (2.55) 

For u ^ 1 the integral in eq. (j2.47p is dominated by the region u' ^ u ^ 1 , and we note that 

n[a (u - u')] "'=" a{u- u') . (2.56) 

In this limit we obtain 



6 r njaju-u')] 1 6 _ ^^ 

a Jo [l-W]^ (i u) an ^(^ + i)U u) 



(2.57) 



Inserting this result in eq. (|2.47p we find two solutions: f3 — 2, —3 which are precisely the singular and regular solutions 
found in the long- wavelength limit [l|. Therefore, the dominant behavior for u — > 1 for all values of a is 

s «^i A(a) ,, , a(u) 

<5(fc,u)"=^-^^-A(a)^^ (2.58) 

$(fc,u)"=^ A(a) (2.59) 

where in order to find the function A{a) we have to integrate the Gilbert equation from u = 0, where the initial 
conditions were defined, up to w ^ 1. <&(fc, u) is a slowly varying function of u in the sense that it is always finite as 
we see from eq. (|2.59p . 



The transfer function: T(fc) is defined from the limit [i| 

aeq 6{k,u) 



1^"^ -t\ Jn^ = 1^"^ '^(^' ") ' (2.60) 

ti^i a(u) o(fc,0) M^i 
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where we have used the relation eq. (|2.46p and the initial condition eq. (|2.48p . It is customary to normalize the transfer 
function so that T{k) ^ 1 as fc ^ 0, obtaining 

r(fc) ^ hm |i^ = I hm (1 - uY Sik, u)^l $(fc, ^=.1)^1 A{a) . (2.61) 

A numerical study of (5(fc,u) and T(fc) depicted in figs. [TII51 [71 [MTPl show that T(k) decreases with fc (or a). This is to 
be expected, k — corresponds to the mode that grows the fastest as it is the most gravitationally unstable, larger 
values of k result in slower growth as a consequence of free streaming. The behaviour of the k = mode is that of 
cold dark matter eq. (|2.52p . 

Therefore, T(k) is a measure of the suppression of the k > wave modes and the fc-scale of suppression of T(fc) is 
necessarily related to the free-streaming wavenumber. 

The final power spectrum Pf{k) is related to the initial one Pi{k) as [l| 

Pfik) = THk)P,ik). (2.62) 

If perturbations do not grow or decay substantially during the prior, radiation dominated phase, Pi{k) is nearly the 
inflationary primordial power spectrum [l|. 

Gilbert's equation (|2.47p can be solved in a power series in u, appropriate for short times. We find from cq. (|2.43p 

6^d,o Jo P^ dpFi{k,p;0) 
Thus, we find from eqs.dMSD and (PT7)) . 

S{k, u) "=" 1 + u^ [3-ba^]+0 (u^) (2.64) 

Therefore S{k, u) starts out by growing or decreasing for small u depending on whether a < Uc or a > ac, respectively, 
where 

3 

Namely, for a < ac the mode 6{k, u) is gravitationally unstable from the start {u — 0) while for a > ac the mode 
starts being gravitationally stable. However, due to the cosmological expansion all modes redshift and eventually 
become gravitationally unstable at some time u < I. As can be seen from ea. (|2.58p . for u ^ 1 all modes are unstable. 

At early times, ac determines the stability limit of density perturbations. 

We carried out an exhaustive numerical study for cold, hot and warm dark matter (Maxwellian (MB) , BE and 
FD distributions). We find for large a that T{k) decreases exponentially in the Maxwell-Boltzmann case while for 
both bosons and fermions T{k) asymptotically decreases with a power law or exponentially (see below sees. IIV Al and 
HVBJ) . 

At different scales, the suppression of T(k) is described by a characteristic wave vector kchar which depends on the 
initial conditions, the particle statistics and the regime of wavenumbers. 

We obtain kchar from the numerical study for different cases and k regimes, in all cases considered we find that 

^char ^^ '^ f sy^eq) • 

C. Thermal and Gilbert's initial conditions 

It remains to specify the initial perturbation Fi{k,p;u = 0) of the distribution function. Although in principle 
this function should be obtained from the full evolution through the prior radiation dominated stage, in what follows 
we consider the physically motivated case of adiabatic perturbations for which the perturbation Fi(k,p]u = 0) 
corresponds to a temperature perturbation, 



namely, 



Td,0 



p., n^ rr. dfo{p,T) ATjk) 

Fi(fc,p;u = 0) =T — — j^ . (2.65) 
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These are the initial conditions proposed in refs.[2l|, l23|- Alternative initial conditions were proposed by Gilbert [2g |. 
who chose 

Fi(fc,p;u = 0) = /o(p)C(fc), (2.66) 

with some unspecified fmiction C{k). 

D. Solution of the Gilbert equation in powers of k. 

We can solve the Gilbert equation (|2.47p in powers of a^, [that is, powers of /c^, see eq. (|2.36|) ] 

5{k,u) = &i:>{u) + a^ 5i{u) + 0{a^) , (2.67) 

where So{u) is given by eq. (|2.52p and Si{0) — 6'i{0) = 0. Inserting eq. (|2.67p in eq. (|2.47|) . using eq. p.63p and deriving 
twice with respect to u leads to the equation 

where b = b for Gilbert initial conditions and & = | 6 for temperature initial conditions. This inhomogeneous 
differential equation for Si{u) can be easily solved using the Green's function of the differential operator in the 1. h. 
s. We find for u ^ 1 after calculation, 

Q U 1 

-:r + C(l) Gilbert initial conditions, 

35 (1 — uj'' 

Siiu)''^' { (2.68) 

-7T + 0(1) Temperature initial conditions. 

105 (1 — u)^ 

Using eq. (|2.6ip we now obtain T{k) in powers of a^, 

8 h 
1 — — a^ + 0{a'^) Gilbert initial conditions, 

T{k) = { (2.69) 

OIL 

1 — - a ~\- 0{a ) Temperature initial conditions, 

Do 

where b is given by eqn. (|2.63p . 

E. Gilbert's equation as an integro-differential equation. 

The results obtained above suggest to convert the Volterra equation into an integro-differential equation. Taking 
the second derivative with respect to u of eq. (|2.47p yields 

r,, ^ 6S(k,u) „ 9 /"" , / /"^ , A 7! , ^ sm\a v (u — u' )] S(k,u') ••., , .„ „„n 

S{k,u)- ^ ,2 +6«' / du' dyy^foiy) ^ ^^ ^-^ -^-^^ ^ I[k,u] . (2.70) 

u — uj Jo Jo '^ y u "" J 

It proves convenient to write the following identity in the y-integral: 

y'^y^¥ + yHy^-¥), (2.7i) 

where 



2/2=/ dyy^foiy) with / dy y' My) = 1 , (2.72) 

and to introduce 



o 2 _ 2 — 2 fc Td,o / y2 



m Ho y 3 ilM fleq 



(2.73) 
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where ea. (l2.36p was used. 

In terms of kfs [see eas. ipTTS)) . 1^^ and ([27721)] we find that 



2fc2 



2fc2 



kfgV'eq) ^fs\ ' ^^1 



, kf siO) 



/ 3 JIm „ _]]]_ 
2 7 ° Td,o 



(2.74) 



Using the original eq. (|2.47p to replace the term with y^ in terms of 5{k,u) and I[au\, we obtain the following 
integro-differential equation for density perturbations 



S(k, u) ~ ^P^ + 3 7' S{k, u) - f du' K{u - u') -^^^ - i[k, u] + 3 7' I[a u] 
where the non-local kernel K{u — u') is given by 



(2.75) 



and we note that 



K{u ~u')^Qa j y {y'^ - y^) /o(y) sin[a y {u~ u')] dy 
Jo 



\ m / 



is the three dimensional velocity dispersion of the non-relativistic particles today. 
Using ea. (|2.26p . 7 [eq. (|2.73p ] can also be written as 

V2k 



7 



kf.iteq) x/3 



' U„,„;0)-^""<"^»> 



V^ A 



where A is the wavelength of the perturbation. 

Using eas. (|2?24| . ((03)) and (|237ll . kfs{teq) and lfs{v;0) are given by 



.,.(,.,) ^ 2:5^ (f )* i^ [kpcl- 



A/s(ie9) = a/:7 7i"Z/^(?;;0) = 616 J y2 ( — ) ' kpc 



(2.76) 



(2.77) 



(2.78) 



(2.79) 



(2.80) 



The transfer function for the three different cases (MB, FD and BE) turns to be parametrized by the dimensionless 
ratio 7 = \/2 k/kfs{teq) allowing us to compare the three different cases, in fig. [TOl below. 

The alternative form eq. (|2.75p of the original Volterra (Boltzmann-Vlasov) [eq. (|2.47p ] equation has several merits : 

• Neglecting the non-local term and the inhomogeneity, eq. (j2.75p can be written in a more familiar form in terms 
of cosmic time t: 



d^S ^^dS 



a\t) 



AttG pM{t) 



(5 = 



(2.81) 



where {V'^) is given by eq. (|2.77p and {V'^)/a'^{t) is the physical squared velocity of the non-relativistic particles. 
The term between brackets in ea. (|2.8ip can be written as 



a^t) 



k' k}s{t) 
a^{t) a^{t) 



(2.82) 



where kfs{t) is the comoving free streaming wavevector ea. (|2.22p . This equation must be compared to the usual 
fluid equation [24] in which the comoving Jeans wave vector kj{t) replaces kfs{t) with 






(2.83) 
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where Cg is the comoving speed of sound during matter domination. This expression shows that the comoving 
Jean's wavevector in the fluid description scales as \/a(i) during matter domination, just as the free streaming 
wavevector fc/s(i) does. Thus, we see that the long- wavelength limit a — > (7 ^ 0) of eq. (|2.75|l yields the 
familiar fluid description. 

Therefore the contribution from the non-local kernel K{u — u') is a correction to the fluid description. It will 
be seen below that this correction becomes important at small scales k > kfs{teq)- 

• In the long wavelength limit a — > (7 ^ 0) the kernel K{u — u') decreases as 7"^ <C 7^ and its contribution to 
the dynamics of density perturbations becomes negligible. 

Furthermore, expanding the kernel K{u — u') in a power series in k it can be seen that each term corresponds 
to higher correlations of the form 



©'^") -(©'■)(©') 



These correlations are typically neglected in the moment-hierarchy of the Boltzmann equation [l| . 

• At short times u ~ it follows that 


hence, its contribution is subleading for u < I/7 since (5(fc, m ^ 0) ^ 1 by the initial condition. 

• As w ^ 1, the u' integral in the kernel is dominated by the region u ^ u' ^ \ for which 5{k, u') ^ A(a) j (1 — u'Y 
as obtained in ea. (|2.58p . therefore it follows that for u ^^ 1, 

^ du' K{u ~ u') i^^ "=^ a' [7 - i¥f] In ^ « 7-^ . (2.84) 

(1 — u')^ 1 — u (1 — U)'' 

Therefore in this region the contribution of the kernel is subleading as compared to the first three terms in 
eq. p.75|) although it could be larger than the inhomogcneity on the right hand side. 

• It explicitly yields the asymptotic behavior S{k,u — > 1) ex 1/(1 — u)^ as a consequence of the second term in 
ea. (l2.75p and the analysis above. This feature will be seen manifestly in the analysis below. 

F. The exact transfer function and its Fredholm expansion. 

The above analysis points out that in all the regimes of interest the non-local contribution from the kernel is small as 
compared to the first three terms in the eq. (j2.75p and can be treated perturbatively. For this purpose, it is convenient 
to write ea. (|2.75p as 



6 

(1-u) 

The source term is 



^('=' ") - 77^2 '^(^' w) + 3 7^ S{k, u) = S[6; u] . (2.85) 



S[S;u] = Sb[u] + Snb[S;u] (2.86) 
where S'b[u] and Snb[S]u\ play the role of the Born and next to Born approximations, 

Sb[u] - 1 + 372/ (2.87) 

Snb [S;u] = r du' K{u-u') .f^'"2 ■ (2-88) 

As discussed above [see ea. (|2.81|) ]. in cosmic time the left hand side of eq. (|2.85p is precisely of the form of the Jeans 
equation for fluids. The right hand side iS'[(5;u] may be interpreted as an additional pressure term non-local in time 
which includes free-streaming through the initial condition in Sb- In the long wavelength limit S'[(5;'u] vanishes 
along with the usual pressure term 87^ 6(k,u) leading to the usual fluid-like equation for cold-dark matter density 
perturbations. 
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Eq. (|2.85p lends itself to a Fredholm (iterative) solution, the basis of which is the homogeneous solution and the 
retarded Green's function of the differential operator on the left in eq. (|2.85p . 
For S'[(5, m] = 0, eq. (|2.85p becomes the homogeneous differential equation 



P'et, 



6 



S'^"\k,u)+3j^ 6^"\k,u) = 0. 



(1 - u)2 
which is solved in terms of Bessel's functions. The general solution is given by 

,5(o)(z)=A/ii(z) + B/i2(^), 
with 

z = zq (l-u) , zo = V37 , 
hi{z) and h2{z) are related to the spherical Bessel functions ^2(2:) and J2{z) [49|, 

3 . 

cos z H — sm z 



hi{z) = ~z niiz) = ^ - 1 



h2{z) = z J2iz) 



3 \ . 3 

— — 1 I sm z cos z 



(2.89) 

(2.90) 
(2.91) 

(2.92) 
(2.93) 



hi{z) and h2{z) are the fundamental irregular and regular solutions, respectively, with the small argument behavior 
for z ^ (m ^ 1): 



hi{z) 



Z-+0 



"2 U = 



and Wronskian 



W\z\ = /i^(z) /ii(z) - /i'i(z) h2{z) = 1 , 



(2.94) 



(2.95) 



where the prime denotes d/dz. The mode functions /ii,2(w) are identified as the growing and decaying solutions of 
Jeans's fluid equations (|2.8ip . The initial conditions ea. (|2.48|) yield 



Using eQ. (|2.95p we find 



A^h'2{zo) , B = -h\{za) 



(2.96) 
(2.97) 



and the homogeneous solution is given by 

<5(")(z) = /i^(zo) hi{z) - h[{zo) h2{z) . (2.98) 

Using the small argument behavior eq. (l2.94p . for 7 = (fc = 0) the homogeneous solution 5^^\z) reduces to eq. p.52p . 
The general solution: The inhomogeneous equation ()2.85p can be formally solved in terms of the retarded Green's 



function of the differential operator on its left hand side which obeys 

^ d^ 6 



+ 37" 



G[u,u') = 5{u-u') 



dv? (1 — uY 

G(u,u') can be explicitly written in terms of hi_2{u) [eqs. (|2.92p - (|2.93p ] as, 

1 



G(u,m') 



V3 



7 



[hi{u) h2{u') - h2{u) hi{u')] Q{u - u') 



The formal solution of eq. (|2.85p is then given by 

1 



(5(fc,u) = (5(°)(z) + 



V3 



7 "'o 



[hi{u) h2{u') - h2{u) hi{u')] S[5\u'] du' . 



(2.99) 



(2.100) 



(2.101) 
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It is convenient to separate explicitly the source term Sb[u] ea. (|2.87p which does not involve the kernel K{u — u'), but 
only involves the free streaming solution of the Boltzmann equation in absence of self-gravity. Integrating by parts 
twice the term with /, using the differential equation (|2.99p obeyed by the Green's function and the initial conditions 
ea. (|2.44p . we obtain 

S{k,u) = (5(i)(fc,u) + -^^ / [hi{u) h2{u') - h2{u) hi{u')] Snb[5;u'] du' , (2.102) 

V 3 7 Jo 

where 

(5(i)(fc,u) - I[au] + -^ r [hiiu) h^iu') - h2{u) hi{u')] ,/^""'j„ du' . (2.103) 

V37J0 {l~u'Y 

The gravitational potential <^(k,u) can be analogously expressed using its relation with the density perturbation 5 
eq.(I23Sl) and eq. (|2.10ip with the result 

^k,u) = ^'^^\k,u) + -^{l-uf f [hi{u)h2iu')~h2iu)hi{u')] SNB[S\u']du' , (2.104) 

V3 7 Jo 

where 

$(i)(fc,u) = (l-u)2,5(i)(fc,u) , (2.105) 

is the Born term. 

It is straightforward to check that for fc = the second term in ea. (|2.104p vanishes and $'^^(0, u) is given by 
eg. d^TMl) (for this note that /[O] = 1). 

Free streaming for fc 7^ leads to collisionless Landau damping of the gravitational potential and density perturba- 
tion. Although density perturbations still grow as 1/(1 — w)^ for u -^ 1 the gravitational potential is bound |$| < 1 
and is a slow variable. 

We obtain an exact expression for T(k) inserting eqs. (j2.104p into the definition of the normalized transfer function 
ea. (|2.60p and using the small argument limits eq. (|2.94p . 



^''' = ^i''*> 



I[au] ,1 

onb[o;u\ 



(l-u)2 6 



du , (2.106) 



where 5{k,u), argument of Snb [see ea. (|2.88p ]. is the solution of the integral equation (|2.102p and we have neglected 
terms proportional to Ueq ^ 10^^. In obtaining this result we have used the fact that I[au] given by ea. (|2.43p is finite 
at M = 1 and that in this limit 5(k,u) oc 1/(1 — u)^. Only the terms featuring hi{u) outside the integrals survive in 
this limit. The terms with hi(u') inside the integrals yield at most terms ex 1/(1 — uY but they are multiplied by 
(1 — w)^, and vanish for u — > 1. Furthermore, in the long wavelength limit J[0] = 1 ; Smb oc 7* and from the small 
argument behavior eq. (l2.94p it is straightforward to confirm that T{k = 0) = 1. 

The integral equations (|2.102p and (|2.104p can be iterated generating the usual Fredholm series as 

6{k,u) = 6^^Hk,u) + 5''^\k,u) + --- , 

$(fc, u) = $^^' (fc, u) + $(2^ {k, u) + • • • , (2.107) 

where 6^^'{k,u) is given by eq. (|2.103p to which we refer as the Born term because of its obvious similarity to the 
solution of a potential scattering problem, 

1 /■" 

5^'^\k,u) = -^ [hi{u)h2{u')~h2{u)hi{u')] SNB[S^^'^]u']du' (2.108) 

V 3 7 "'0 

Snb[6^^^;u] = r du' K{u~u') '^''^^''^2 , (2.109) 

$(^)(fc,M) is given by eq. (|2.105p and 



$(2)(fc^u) = {l-uf -^- \ \hi(u) h2{u')~h2{u) hi{u')] Snb[6^^'^;u'] du' . (2.110) 

V 3 7 Jo 
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The small argument behavior eq. (|2.94p makes manifest that the solution S{k,u) indeed has the asymptotic behavior 
oc 1/(1 — u)^ as u — *■ 1. Furthermore, since the function /[au] describes the free-streaming solution of the Boltzmann 
equation in absence of self-gravity, the formal solution eq. (|2.102|) exhibits explicitly the free-streaming phenomenon. 



The transfer function has also a systematic Fredholm expansion which up to second order yields 



r(fc) = 



10 



V3 



7^ "'0 



/i2(u) du 



(T^^ + 6^^^['^ 



TB{k)+TNB{k) , 



where we refer to the first order term 



TBik) 



10 



\/3 T' Jo 



n2(u) du 



(1-U)2 ' 

as the Born term because its origin is the Born term for the gravitational potential ^^^' (fc, u) and 



TNB{k) = 



3^3 



7" 



/i2(w) Snb[5 ';u]du , 



(2.111) 



(2.112) 



(2.113) 



is the next to Born (second order) correction, where (5^^^(fc;M) and S n b[S^'''^ ', u] are given by eas. (|2.103p and (|2.109p . 
respectively. 

At this point it is worthwhile to emphasize the following important aspect: the first order. Born term in the transfer 
function only depends on the initial condition and describes the free-streaming suppression. The second order term, 
the integral of Snb, contains the information on the higher moments of the distribution function of the decoupled 
particles and the corrections to the fluid approximation through the non-local kernel K{u — u') eq. (|2.76|) . The details 
of the distribution function beyond the first moment (p^) enter the transfer function at second order and beyond in 
the Fredholm series. It will be shown below that these contributions include memory of gravitational clustering and 
become important at short wavelengths 7 > 1. 



Eq. (|2.11ip provides a useful approximation to the transfer function for general initial conditions and distribution 
functions. Specific examples for which we obtain an exact numerical evaluation (see below) show that the second 
order approximation ea. (|2.11ip is remarkably accurate in a wide range of scales. These are some of the main results 
of this article. 



III. COLD DARK MATTER: NON-RELATIVISTIC, MAXWELL-BOLTZMANN DISTRIBUTION. 



The unperturbed distribution function for particles that decoupled while non-relativistic in local thermal equilibrium 
is a solution of the coUisionlcss Boltzmann-equation (in absence of gravitational perturbations), given by 



f^{Pf{t),t)^Me 



-T- 



iTjc 



(3.1) 



where the explicit expression for the normalization factor M may be found in ref.[39|, Pf{t) is the physical momentum, 
a{t) is the scale factor and T^, a^ are the temperature and scale factor at decoupling respectively. Since p = Pf{t) a{t) 
is the comoving momentum (equal to the physical momentum today) and Td ad = Tdfi is the decoupling temperature 
today, the unperturbed Maxwell-Boltzmann distribution function can be written in terms of dimensionless variables 
as 



/o(2/;a;) =7Ve 2- 
and the normalized distribution function is given by 

Uv) = 

j/2 [eq. (|2.72p ] and a [eq. (|2.73p ] are in this case. 



y 



_P_ 
Tdfi 



m 
T'd 



^[2x\ 



(3.2) 



(3.3) 



y" 



Td 



(3.4) 
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The typical values of the masses and decoupling temperatures for WIMPs are m ~ 100 GeV and Td ^ 10 MeV pa |. 
for which g^ ~ 10 ^2]. Using the result for the free-streaming wavevector today ea. (|2.79p we find 

kfsite,) = ^ (f ) ^ (loO^) ' (loS^) ' • (^-^^ 

Within the cosmologically relevant range of scales where the linearized approximation is valid, eq. (j2.39p and using 
([273)) it follows that in this range 7 < 10^*". 

The kernel n[a {u — u')] that enters in the Volterra equation (j2.47p is best obtained by performing the momentum 
integrals of the distribution function in cartesian coordinates rather than performing the angular integration first. 
This is a consequence of the fact that the distribution function is a function of p ^ . Wc obtain the following equation 
for density perturbations 

5{k,u)~Q r du' (u - u') e-^ <"-"')' .'^^^'"2 ^ -^[7"] ■ (3-6) 

where /[7U] [ea. (|2.43p ] is the inhomogeneity determined by the initial condition. Taking two derivatives with respect 
to u and using eq. (|3.6p we obtain 

5{k,u) ^5(fc,u)+3 72,5(fc,u)-6 74 Tdu' (u - u')^ e^^ <"-"')' ^^^^^ = /[7 u] + 3 7^ /h u] . (3.7) 

1 - u^ Jq [1 - u'Y 

From this expression we obtain the explicit form of the kernel Kiu — u') that defines the source term ea. (|2.88p and 
enters in the transfer function eq. (|2.106p . namely 

K{u - u') = 6 7* (u - u'f e"^ ("""')' . (3.8) 

In the cosmologically relevant case 7 < 10~^, we can safely approximate the density perturbation, gravitational 
potential and transfer function by their Born term in the Fredholm series, namely 

<5(i)(fc,u)=/[7 «] + -£- r [hiiu) h2iu') - h2iu) hi{u')] jT-^du' , (3.9) 

V 3 7 Jo U " "" J 

$(i)(fc,u) = (l-u)2,5(i)(fc,u), (3.10) 

Tsik) = -^ f h^iu) -M^ du . (3.11) 

For the sake of comparison with the exact result and to display the dependence on initial conditions we study both 
initial conditions eas. (|2.65p and (|2.66p . The calculation of I[ju] in both cases is best performed with cartesian 
coordinates, we find for the case of temperature perturbations ea. (|2.65p 



1 1 2 2 
I--7 u 



(3.12) 



It [7 ""] 

and for Gilbert's initial condition eq. ()2.66p 

lG{-fu]^e-i'''''\ (3.13) 

They are connected by the simple relationship 

lTbu]=(^l + r^j-^^lGbu]. (3.14) 

Figs. [T] and [H display the exact results for ^{k;u) and T'^{k) obtained by numerical integration of eq. (|3.6p and 
using the relation (|2.46p . compared to the solution obtained with the Born approximation. It is clear that the Born 
approximation is remarkably accurate for 7^1, namely k < kfs{teq)- Fig. [3]compares T'^{k) for the initial conditions 
eq. (|3.13p and eq. (|3.12p . The exact numerical results are indistinguishable from the Born approximation in the range 
displayed. 
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FIG. 1: Non-relativistic, Maxwell-Boltzmann distribution. Left panel: The gravitational potential ${k,u) vs. u for 7 = 
0, 0.25, 0.4, 0.75, 1. Right panel: the transfer function T^{k) vs. 7. In both cases the solid line corresponds to the exact 
solution of eq. (|2.47|) with the initial condition ea. (|3.13p and the dashed line to the Born approximation eas. (l3.10|l and (|3.11[l . 
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FIG. 2: Non-relativistic, Maxwell-Boltzmann distribution. Left panel: <l>(fc, it) vs 7 for 7 = 0,0.25,0.4,0.75,1. Right panel: 
T'^ik) vs. 7. In both cases the solid line corresponds to the exact solution of eq. (I2.47|) with the initial condition ea. (l3.12p and 
the dashed line to the the Born approximation eqs. (|3.10[) and 1)3. 11[) . 



A. Short w^avelength: memory of gravitational clustering 



Although for WIMPs only the region of 7 <C 1 is of cosmological relevance as discussed above, it is important to 
study the opposite limit 7 ^ 1 to understand how the memory kernel K{u — u') modifies the transfer function. It 
proves convenient to change variables to z = 7 (u — u') and to replace in eq. (|2.88p the density perturbation 5 by the 
gravitational potential <i>(/c, u) which is bounded in time, to write 



Snb [u) = 6 



z^e-^ L ^1 



[l--+f]' 



dz 



(3.15) 



Since <I>(fc;u ~ 0) ~ 1 it is clear from this expression that Smb is negligible either at small time (w <C I/7) or in the 

long wavelength limit 7^1. The factor z^ e~~2" in the integrand grows as z^ at small z, attains a maximum at 
z = Vs, namely u~ u' ^ I/7, and falls-off exponentially for z > \/3 for 1 — u 3> I/7. 

Therefore, Smb begins to contribute for u > I/7. Since $(fc,u) is still ~ C(l) during this interval it follows that 
Snb becomes of 0(1), whereas for large 7 the Born term for Gilbert's initial conditions, 



Sb.c 



7 



[2- 



-I'u^] 



7 u 



(3.16) 



is suppressed for u ^ I/7. Analogous conclusions follow for temperature perturbations where the expression for Sb.a 
follows combining eqs. (|3.14p and (|3.16p . 
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FIG. 3: Non-relativistic, Maxwell-Boltzmann distribution. Comparison between T (k) vs. 7 for the initial conditions eas. (|3.12p 
and 1)3. 13p . The exact result and the Born approximation are indistinguishable in this range. 



Therefore, for 7^1 there is a crossover in behavior: during < u < I/7 the Born term which dominates is of 0(1) 
and the second order correction is negligible, while at m ~ I/7 both terms become of the same order, and for u > I/7 
the second order correction becomes important. Finally, for 1 — u <C I/7 the small z region dominates yielding a 
logarithmic behavior for the second order correction. When this correction becomes important for 7^1, the transfer 
function T{k) is very small, suppressed at least by the prefactor 1/7"^. 

For 7 3> 1, in the region with 71* 3> 1 the upper limit of the integral can be taken to infinity. For I/7 <C w <C 1 — I/7 
the dominant contribution arises from the region u ^ where $ '--^ 1 is largest, in this region the full integral is of 
0{1). 

Furthermore, since to leading order $ decays via free-streaming on a time scale u ^ I/7, its derivative is large 
during the initial free-streaming regime but becomes small for u > I/7 as can be gleaned from Figs[T][2l This analysis 
leads to the conclusion that S^Biu) becomes important for 7 3> 1 in the region 7 m 3> 1, where it is of the same 
order (or larger) than the free streaming contribution (Born term) . When Snb (u) becomes non- negligible the transfer 
function is suppressed by ^ 1/7'^. Therefore, the corrections to the Born term become relevant at small scales when 
the transfer function has diminished substantially. 

In the region u 3> I/7 the following Markovian approximation is reliable: 



$ [fc; M ] « <i> [fc; m] - $ [fc; m] - 



(3.17) 



We emphasize that this Markovian approximation is valid only after the initial free streaming transient for u > I/7, 
since during the transient $ '^ 7. Therefore, although such Markovian approximation is available after the initial 
transient, the value of the gravitational potential must be matched to that at the end of the free-streaming period in 
order to provide the full dynamical evolution. 

The analysis above suggests that the influence of the memory kernel K{u — u') becomes important for short 
wavelengths 7 > 1. Thus, for scales A < Xfs{teq) ^ lfs{v', 0) where 7 > 1 we need to include the second order term in 
the Fredholm expansion in the transfer function given by eq. (|2.11ip . where 



Snb[S'^' 



6 7^ / du' {u - u'f 



(„_,). 6('Hk,u') 



[1 - u'] 



'12 



(3.18) 



and S^^'{k,u) is given by ea. p.9p . Higher order terms are further suppressed by extra powers of I/7 in the transfer 
function. 

The addition of the second order correction yields a remarkably accurate fit to the exact solution of the Boltzmann- 
Vlasov equation in a wide range of scales down to scales A <C \fs{teq)- 

Fig. U displays the comparison between the exact result for T'^{k) compared with both the Born approximation 
and the next to Born approximation for initial conditions eq. (|3.13p . Fig. [5] displays the same comparison for the 
initial condition eq. ()3.12|) . These figures confirm the analysis presented above: the Born term is a fairly accurate 
approximation in the region 7 ;^ 1 and the second order correction becomes significant at 7 ~ 1 . Including the second 
order correction yields a remarkably accurate approximation in a wide range of scales < 7 < 5. 
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FIG. 4: Non-relativistic, Maxwell-Boltzmann distribution. Left panel: The transfer function T (k) vs. 7 for the exact solution 
and the Born approximation. Right panel: T'^{k) vs. 7 for the exact solution and the Born approximation plus second order. 
Gilbert's initial condition ea. (l3.13p was used. 
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FIG. 5: Non-relativistic, Maxwell-Boltzmann distribution. Left panel: the transfer function T^{k) vs 7 for the exact solution 
and the Born approximation. Right panel: T'^(k) vs. 7 for the exact solution and the Born approximation plus second order. 
Initial condition for temperature fluctuations eq. (|2.65p 



We obtain T(fc) for large scales expanding eq. p.lSp for small 7U 

r ,7 u^O 1 2 2 , /o/ 4 4\ 

iGnul = I--7 w +C(7 Mj 



(3.19) 



and follow the same steps leading to eq. (|2.69|) . setting 6=1/2 according to eqs. (l2.63p and (I3.19p . and replacing a by 
7 in ea. (|2.69p . with the result 



where 



T(fc) = 1 - 



70 



Oh') 



(3.20) 



70 



f = 1.62. 
/^ = 1.42. 



Gilbert initial conditions, 
Temperature initial conditions. 



(3.21) 



70 characterizes the fall-off of T{k). We find for a given initial condition 70 is independent of the statistics for thermal 
relics. 

For WIMPs the long-wavelength approximation (|3.20p describes T(k) in the whole range of scales of cosmological 
relevance for structure formation in the linear regime, since for these scales 7 < 10^^ and the Born approximation 
describes the DM density and gravitational perturbations outstandingly well. 
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Although for the scales of interest for structure formation T{k) given by (|3.20p gives the correct description, it is 
illuminating to study the small scale behavior of T{k). 

The precise numerical solution of eq. (|3.6|l shows that T{k) for temperature initial conditions has a zero at 7 = 7.7 . . . 
as shown in fig. [TUland is negative for 7 > 7.7 . . . decreasing exponentially for 7 > 10 as 

|r(fc)| -e^V^^S^J (3.22) 

where xmb — 1.3 and ^mb — 1-55 both for temperature and Gilbert initial conditions cqs. (|2.65p and (I2.66p . 
We draw the following important lessons from this study: 

• For k < kfs{teq), the Born approximation gives a simple and very accurate description of the transfer function. 
It inputs the free streaming solution in absence of self-gravity and the mode functions associated with the fluid 
description. These only input the average squared velocity, therefore the first moment (p^) of the distribution 
function. This approximate but fairly accurate solution yields a simple tool to understand different distribution 
functions and initial conditions in generic situations. 

• The short wavelength region of the transfer function k > kf^iteq) is remarkably well described by correcting 
the Born approximation with the second order term in the Fredholm solution. This correction incorporates 
higher moments of the distribution function and includes important memory effects of gravitational clustering 
and corrections to the fluid description. This approximation while just slightly more involved than the Born 
term also yields a rather simple and systematic tool to study the effects of higher order correlations and the full 
structure of the distribution function along with memory of gravitational clustering in generic situations. 

• We see that the scale characterizing the suppression of T{k) with 7 decreases with k (and 7). Namely, 70 
characterizing the suppression for small k is larger than the characteristic suppression scale for larger k. 

Nevertheless, the numerical study shows that the characteristic suppression scale kchar related to a characteristic 
dimensionless ratio "fchar by eqs. ()2.33p . (|2.36|) . ()2.38|) and (|3.4p . namely: 



Knar = \ Ho v/^W^ (f ) ' -^ 7cw = 4.17 (I) ' (j^) ' {-^) ' ^ . (3.23) 

is such that in all cases considered we find jchar ^ C'(l). Therefore, kfs{teq) obtained from free particle 
propagation [eas. (|2.26p - (|2.27p ] and given by equations (|2.79I3.5P only differs by a factor of order one from kchar 
[eq.JSSl]. 

IV. WARM AND HOT DARK MATTER: FERMIONS VS. BOSONS 

In this section we consider thermal relics that decoupled in LTE while ultrarelativistic, but that have become 
non-rclativistic during matter domination. Their normalized distribution function after freeze-out are given by 

Uy) - ^ ^ FD , (4.1) 

/o(.)^^^ BE, ,= ^, (4.2) 

for FD and BE respectively, with Td.a being the decoupling temperature today. We study each case separately. 

It is convenient to state the following general result in order to estimate the large momentum contribution of the 
various integrals. Consider a generic integral of the form 

F{y) sm[yz] dy , (4.3) 



where F(0), -F'(O), F"(0) • • • are finite and F{oo) = 0. Integrating by parts consecutively we find the large z behavior 

j^ F[y) sin[yz] dy ^ ^ ~ ^ + O [^^^ ^^.4) 

This result will be useful in the analysis that follows. 
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T(fc) decreases in k with a characteristic scale kchar which depends on the initial conditions, the particle statistics 
and the regime of wavenunibers. kchar translates into a characteristic scale achar in the dimensionless variable a. 
kchar is related to Uchar by eqs. (|2.33p . (|2.36p and (|2.38p : 



kchar — -^ Hq y/^M acq ( "^ ) 



^char — U 



.00417 achar (f)' i^[kpc]-' 



-t C7nb 

For fermionic and bosonic thermal relics achar ~ 0(1) because in these cases y^ ^ 0(1)- 



(4.5) 



A. Fermions 



For fermions that decoupled relativistically in LTE with the normalized distribution function eq. (|4.1|) . it follows 
that 



y 



2 _ 



y^ foiy) dy = 15 



C(5) 



12.939- •• 



/o C(3) 

leading from cq. (|2.79p to the free-streaming wavevector today 

fc/,(t,J =0.00284... (I) ^j^[kpc]-i 



For the initial condition eq. (|2.65p the inhomogeneity is 



It [a u] 



y^ e^ sm[y a u] 



dy 



9 C(3) Jo (e^ + 1)^ c^u 
whereas for the initial condition eq. (|2.66p 

2 f°^ y sm[y a u] 



E 



\n+l 



1-1 



3 C(3) ;^^ (n2 + z2)2 [ 3 (n^ + z^) 



3 C(3) Jo ey + 1 au 



^^ 3C(3)5^(n2 + z2)2 



z = a u 



(4.6) 



(4.7) 



(4.8) 



(4.9) 



Using the result eq. (|4.4p we find that for both initial conditions the asymptotic behavior of the inhomogeneity for 
a M ^ 1 is given by 



Ig [a u] 



a. u— >oo 



3 C(3) (a m)4 



O 



1 



, It [a u] 



a u — »oo 



1 



9 C(3) (a uY 



o\ 



1 



(4.10) 



which must be contrasted to the Maxwell-Boltzmann case for which /^[au] and It[(xu] eqs. (l3.12p and (|3.13p decay 
exponentially. 

Therefore, FD statistics results in free-streaming solutions in absence of self-gravity /^[au] and /T[aM] that fall 
off much slower, with a power law in time. Such long-range feature is also present in the kernel K{u — u'). This is 
unlike the Maxwell-Boltzmann case where such free-streaming solutions eqs. (|3.12p and p.l3p fall-off exponentially on 
a time scale I/7. 

Implementing the result eq. (j4.4p we find that the kernel K{u — u') eq. (|2.76p with the normalized FD distribution 
function eq. (|4.ip falls off as 



I a («-u')-+oo 30 C(5) 

K[u — u ) — 



eCi) [a{u-u')Y 



FD 



(4.11) 



Fig. [5] displays both A'[z]/[6 a] and z^ K[z\/[& a]. Again this case must be contrasted with the Maxwell-Boltzmann 
case eq. (l3.8p which decays exponentially. 

Thus, whereas for 7^1 the Maxwell-Boltzmann distribution leads to a short range memory kernel, the FD 
distribution yields a long range kernel that keeps memory of the initial state and the initial value of the gravitational 
and density perturbations. 

For the FD case eqs. (P?75)) and (|i^ yield 



K(3) 
5C(5) 



7 = 0.482... 7 



(4.12) 
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z = a (u-u') 

FIG. 6: The kernel K{u — u')/[6 a] [see ea. (|2.76p ] vs. z = a {u — u) for Fermionic thermal relics. 



We can now study the transfer function as a function of 7 comparing the exact result with the Born approximation 
ea. (|2.112p and the next to Born term eq. (|2.113p . The analysis presented in the previous section of the different 
contributions also applies to this case: 

For long wavelengths 7 < 1 the second order correction 5[5^,u] is much smaller than the Born term because the 
kernel K{u — u') is of order 7^, hence T{k) is dominated by the Born term and by the free streaming solution. 

At shorter wavelengths 7 > 1 the kernel K{u — u') begins to contribute at u > I/7 and peaks at u' ~ m as indicated 
by Fig. inibut by this time both the density perturbation and gravitational potential inside the integrand have decayed 
substantially leading to a suppressed contribution to the transfer function. 

This analysis leads to conclude that for 7 < 1, T(fc) is dominated by the Born term while second (and higher) 
order corrections become relevant for 7 > 1. 

We have confirmed this analysis by comparing the exact numerical solution of Gilbert's equation with the Born 
approximation and the second order correction to the transfer function for both initial conditions eqs. (|4.8p and (j4.9p . 
The results are qualitatively the same in both cases and we present them for temperature initial conditions in fig. 
[71 The agreement between the exact result and the Born plus second order correction is remarkable. The numerical 
study confirms the analysis above, the second order contribution, which includes the memory of gravitational clustering 
becomes important for 7 > 1, namely at small scales, but when it begins to be appreciable, T{k) has been highly 
suppressed. For 7 > 5 wc find that T'^{k) < 10^^. 

T{k) for large scales follows from ea. (|2.69p and 

for fermions. We find that the large scale approximation for T{k) ea. (|3.20p is also valid in the FD case. 

The precise numerical resolution of eq. (|3.6p for fermions shows that T{k) for Gilbert initial conditions [eq. (|2.66p ] 
and 7 > 4 decreases as 

T(fc) 0, f ^"j '' , xfg :^ 7.6 , 7/s ^ 3.7 , (4.13) 



7 / 



see fig. [TOl 

The behavior of T(fc) for temperature initial conditions [eq. p.65p ] is more involved as displayed in fig. 1101 we find 
that T(fc) decreases exponentially for 4 < 7 < 15 as 



T{k)~cfae -/" 



Cfa ~ 7.6 , 7/al ^ 0.89 . 



(4.14) 



For 7 > 15 T{k) decreases faster than in eq. (|4.14p . it vanishes and becomes negative at 7 ~ 16.72. For 7 > 25 T{k) 
decreases in absolute value as power law: 



T{k) 



2fa2 

1 



, Xfa ^ 5.3 , 7/o2 ^ 0.96 . 



(4.15) 
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We find that the scale of suppression itself slides with scale. 

For Gilbert initial conditions, we see that the characteristic fall-off scale of T{k) increases with increasing k, from 
7o ~ 1.62 ... for small k to -ffg ~ 3.7 for 7 > 6. Instead, for temperature initial conditions, the characteristic fall-off 
scale decreases with increasing k reaching values 7/a — 0.89 — 0.96 for 7 > 4. However, all of these are of the same 
order ^ 0(1) which is a manifestation of a unique characteristic scale ^char ^ C'(l) as also found in the case of the 
Maxwell-Boltzmann distribution. 

Therefore, in terms of wavevector k this observation translates into the statement that the relevant scale for 
suppression of T{k) is kfs{ti,q), although the functional form of T(fc) itself depends on scale and initial condition, in 
this case varying from exponential to power law. 



We anticipate that the power law behavior of the kernel K{u — u') for fermions [eq. (|4.1ip ] leads to a longer memory 
on the initial conditions than in the Maxwell-Boltzmann case eq. (|3.8|) . For 7^1 the range of the Born term and 
that of the kernel are much longer for fermions than for the Maxwell-Boltzmann case. Then, both the free streaming 
solution and the gravitational perturbation (or alternatively the density perturbation) for small values of u' inside 
the integrand in K{u — u') yield larger contributions for FD than for Maxwell-Boltzmann. We indeed find that |r(A:)| 
for thermal Fermions and for a given value of 7 is enhanced relative to that of the Maxwell-Boltzmann particles [see 
fig. [TO]. 




Thermal Fermions Ix(k;u) 
Exact 

Born + 2"" order 




FIG. 7; Fermionic thermal relics. Left panel: T^{k) vs 7 for the exact solution and the Born approximation. Right panel: 
T^{k) vs. 7 for the exact solution and the Born approximation plus second order. Initial condition given by ea. (|4.8p . Similar 
agreement is found with the initial Gilbert's conditions eq. (|4.9|l . 



B. Bosons 



For thermal bosons that decoupled while relativistic with the normalized distribution function ea. (|4.2p we find 



y2 = 12 ^ = 10.352- •• 



leading from eq. (|2.79p to the free streaming wave vector today 

fc;,(ie,) = 0.00317... (|i)^^[kpc]-^ 



For the initial condition eq. (|2.65p the inhomogeneity is 



It [oi u] 



1 



°° iP- pV 



y e^ sa\[y a u] 



dy 



1 



6 C(3) Jo (e^ + 1)^ au 
and for the initial condition eq. (|2.66p 

1 f°° y sm[y a u 



E 



C(3) f^^ (^2 + z2)2 L 3 (n2 + 7J) 



4 z' 



lG[au] 



2 C(3) Jq ev + 1 au 



- 00 



C(3) ^^ (n2 + Z2)2 



z — a u 



(4.16) 



(4.17) 



z = au (4.18) 



(4.19) 



26 

Using the result eq. (|4.4p we find that the asymptotic behavior of the inhomogeneity for a u ^ 1 for both initial 
conditions is given by 



Ig [a u] 



a U— >oo 



2 C(3) {a uf 



o\ 



, It [a u] 



a. u—>-oo 



6 C(3) [a uf 



o\ 



(4.20) 



This is an even slower power law fall off than in the FD case and even much slower than the Gaussian-exponential 
fall off in the Maxwell-Boltzmann case. 

Since /[au] is the free streaming solution in absence of self-gravity (normalized to one at the initial time), we see 
that with either initial condition the Bose-Einstein distribution function leads to a much less efficient free-streaming 
smoothing of the initial perturbation. 

This long range feature associated with the Bose-Einstein distribution is also manifest in the kernel K{u — u'). 
Again, using eq. (l4.4|) we find that the kernel K{u — u') eq. ()2.76|) with the normalized Bose-Einstein distribution 
eq.diSI) falls off as 



, a (u-u')^oo 36 C(5) 

K[u— u) = 



1 



C2(3) u-u' 



BE 



(4.21) 



Fig. [5] displays both iir[z]/(6 a) and z K[z\/{Q a). This fall-off is even slower than in the FD and certainly much 
slower than the Maxwell-Boltzmann case eq. p.Sp . with important consequences explored below. 



zK[z]/6a 




z=a(u-u') 



FIG. 8: The kernel K{u ~ u')/[(3 a] (full line) [see ea. (l2.76p ] and z K[z\/{(i a) (dashed line) vs. z = a (u ~ u) for Bosonic 
thermal relics. 



BE statistics, more specifically the behavior of the distribution function at small momentum, leads to a slow fall 
off with a power law and long range memory, even much longer than the FD case because of the divergence of the 
distribution function as y — > 0. 

The long-range nature of the kernel brings about important consequences: even for u ~ 1 the kernel is sensitive to 
the region u' ^ 0, therefore the initial value of the gravitational potential ^{k, u ~ 0) ^ 1, which is the largest value 
that <I> attains, contributes with a large measure to the integrand. This feature in turn results in that free-streaming 
is less efficient, and a large memory of the initial value of the gravitational potential remains. This produces an 
enhancement in the transfer function as compared to the Maxwell-Boltzmann and Fermi-Dirac cases, in which the 
memory in the kernel is of much shorter range and the initial (maximum) value of the gravitational potential does 
not contribute as much to the integrand for u ^ 1. 

This remarkable difference will be studied explicitly below in a comparison between the three cases. 

For bosons eqs. ((273)) and (|4T6)) yield. 



C(3) 



7 = 0.538... 7, 



(4.22) 



V4C(5) 
and we now have all the ingredients to study the transfer function as a function of 7 to compare to the previous cases. 
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The analysis presented in the previous cases for the magnitude of the contributions from the Born and second (and 
higher) order remains the same. In the long wavelength limit 7^1 the Born term dominates, and the second order 
correction is subleading of order 0{j^). 

For large 7 (small scales) just as in the Maxwell-Boltzmann and Fermi-Dirac statistics, and quite generally, for w ~ 
the integral of the kernel yields a contribution ^ 7^ w^, therefore for large 7 the second order contribution begins to 
be significant at m ~ I/7 but at this time the gravitational potential has decayed significantly for 7 ::^ 1. Although 
the long range memory of the kernel maintains information on the initial values of the gravitational potential, but 
suppressed by a power I/7. Therefore the second order correction becomes significant for large 7 (small scales) but 
is also suppressed by inverse powers of 7, and hence is always perturbatively small. 

It must be noticed that because of its longer range, the second order correction for bosons is comparatively larger 
to that in the Maxwell-Boltzmann and Fermi-Dirac case. 

We have studied numerically the transfer function for both initial conditions and compared the exact numerical 
solution of Gilbert's equation to the Born term and the second order correction. Both cases are qualitatively similar. 
Fig. [9] displays the results of the analysis in the case of the initial conditions corresponding to temperature pertur- 
bations eq. (|4.18p . The results for Gilbert's initial conditions eq. (|4.19p are qualitatively the same, with a remarkable 
agreement between the exact numerical solution and the second order improved Born approximation. 
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FIG. 9: Bosonic thermal relics. Left panel: r^(fc) vs 7 for the exact solution and the Born approximation. Right panel: T^{k) 
vs. 7 for the exact solution and the Born approximation plus second order. Initial condition given by ea. (|4.18p corresponding 
to temperature perturbations. 



T{k) for large scales follows from ea. (|2.69p and 



JG[a^.]"=-^"l-2||(aH)2 + 0([a^,]4) 



for bosons. We find that the large scale expression for T(fc) eq. (|3.20p is also valid in the BE case. 

The numerical solution of ea. (|3.6p for bosons shows that T{k) both for temperature and for Gilbert initial conditions 
[eqs. (l2.65p and (|2.66p ] decreases as a power for 7 > 6: 

, xb — 4.8, "fB — 2.9 for temperature initial conditions , 

xb — 4.6, jB — 3.2 for Gilbert initial conditions . (4.23) 



T{k) ~ ( ^ 



Again the scale of suppression of T{k) increases with increasing k, from 70 ~ 1.62 . . . for small k to jb — 2.9 — 3.2 for 
7 > 6. However, these scales are all of the same order ~ 0(1), evidence of a unique characteristic scale jchar ^ 0(1)- 
Namely, just as in the previous cases, in terms of the wavevector k the relevant scale of suppression is kf^iteq). 

V. SMALL SCALES: STATISTICS AND MEMORY OF GRAVITATIONAL CLUSTERING 



The difference in statistics, namely the distribution function of the decoupled particles, enters in the initial condition 
and in the non-local kernel of Gilbert's equation (|2.47p . The initial conditions are determined from the evolution of 
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perturbations during the cosmological stages prior to matter domination and their dependence on the distribution 
function may be different from the ones studied in the previous sections, while the kernel in ea. (|2.47p is independent 
of the initial conditions. 

The study above revealed a remarkable difference arising from the different distribution functions. We have estab- 
lished that the free streaming solution described by I[a u] has very different asymptotics, falling off exponentially with 
a u in the case of Maxwell-Boltzmann and as power laws in the Fermi-Dirac and Bose-Einstein cases. Furthermore, 
the fall-off is much faster for fermions than for bosons. Since the maximum value of u is u = 1 the asymptotic behavior 
in 7 M actually describes the region 7 S> 1 or scales much smaller than the free streaming scale. 

We have also highlighted that the different statistics lead to an important difference in the contributions from the 
kernel K{u — u'): whereas Maxwell-Boltzmann statistics leads to a short range memory that falls off exponentially in 
{u — u') eg. p.Sp . both Fermi-Dirac and Bose-Einstein lead to long-range memory that falls off with a power law in 
{u — u') [eas. (|4.f fp and (|4.2fp ]. the smallest power, namely the slowest fall-off corresponds to Bose-Einstein statistics 
as a consequence of the infrared enhancement (large population at small momentum) . 



Furthermore, we have shown in the previous sections that the influence of the kernel K{u - 



(the second order 



correction) becomes important at small scales (7 > f ). The long-range of the kernel keeps memory of the early stages 
of the evolution when the gravitational perturbation has its largest amplitude, therefore the long-range nature of the 
kernel leads to an enhancement of the transfer function and the power spectrum at small scales as depicted in fig. If 01 

We can study the small scale behavior in more detail by focusing on the time scales during which the gravitational 
potential varies slowly and a Markovian approximation, such as that discussed in section lTlI Al ls reliable. As discussed 
above and explicitly shown by the numerical study, for 7 3> 1 the gravitational potential decays very rapidly over a 
time scale u ^ f /7, then after evolves slowly until u '^ 1. In order to gain deeper understanding of the dynamics 
during this time scale it is convenient to implement a Markovian approximation directly in Gilbert's equation (I2.47p . 
Upon changing variables in the integrand in eq. (|2.47p a{u — u') = z, ea. (|2.47p becomes 



$(fc,«)-4(l-")' 



n[z] 






dz = {I — u)^ I[au] 



(5.1) 



where the kernel 11 [2;] is given by ea. (|2.4ip . 

For the Maxwell-Boltzmann case [see eq. (|3.6l) ] it is more convenient to change variables to ^{u — u') = z. 

The kernel n[z] vanishes at z = 0, is sharply peaked near z ~ 1 and falls off exponentially for Maxwell-Boltzmann 
(see eg. p. 61) ). as Xjz^ for Fermi-Dirac and as f/z for Bose-Einstein. 

For u 3> 1/a and (I — u) ^ f /a the terms that multiply n[z] may be expanded in a power series expansion in zj a. 
This is the Markovian approximation, in which we obtain 



or du 



<P{k,u) 



n[z] z dz 



i + o{ — 



= {1-uf I[au] 



(5.2) 



The first order equation (|5.2p has a simple exponential solution that determines the decaying behavior of the gravi- 
tational perturbation, where the effective decay rate of gravitational perturbations r[fc;u] is given by 



r[k;u] = I 



6 



a2 (l-w)2 



nfzl dz = I 



fc2 



(5.3) 



The wave-vectors leading to the smaller positive values of the decay rate F are the ones that decay the least and 
for which the power is the largest. This is akin to the criterion that determines the Jeans wave vector in a fluid 
with gravitational perturbations: the Jeans wave vector separates the stable acoustic oscillations from the unstable, 
growing modes corresponding to gravitational collapse. This is also the criterion that determines the free streaming 
wave- vector in Minkowski space-time [37[. 

For Maxwell-Boltzmann and Fermi-Dirac distribution functions the upper limit in the integrals can be safely taken 
to infinity for a u 3> 1, with the result 



n[z] dz ^ [ fo{y) dy = (4 

Jo ^y 



/•OO 

/ z n[z] dz ^ A 
Jo 



(5.4) 
(5.5) 
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where the (finite) constant A depends on the distribution function and the angular brackets stand for the average 
with /q. In these cases we find for K^{t) [see eas. (12.151) and (|2.36p ] 



Klit) 



2 -"OM 



T^ (€;)'(?)= -""^ ""' ^f;)'""' '^■•<»""" 



Remarkably, Kq,{t) coincides with the free streaming wave- vector found in Minkowski space-time 37[, 



(5.6) 



i^*(0) 



4^GpoA/(-^ 



0.563 



V 2 / keV ^ ^ ^ 



(5.7) 



Thus, at small scales the Minkowski result is obtained, consistently with the expectation that at short wavelengths 
an adiabatic approximation to the expansion is reliable [33, |37[ . 

The full transfer function T(fc) cannot be obtained from the Markovian approximation alone and the full study 
presented in the previous sections is necessary. However, it becomes clear that kfs{teq) is not the only relevant scale 
and that there is a further scale K^{tf,q) eq. (|5.7p . The ratio of these two expressions is 



Km 



^^'^ (^> 



y'^ foiy) dy / foiy) dy 



3 Maxwell — Boltzmann 
4.9742 .. . Fermi - Dirac 



Defining the small scale length today 



A$(ie<}) — 



2tt 



K^{teq) 



we find for WIMPs, 



and for thermal fermions, 



A$(teg)~0.5pC ( — )' ( 



100 GeV 



Yi 



lOMeVy 
Td ) 



A$(ie 



939 



kpc (A) 



2 \ 3 keV 



FD 



(5.8) 



(5.9) 



(5.10) 



(5.11) 



The results above are valid for both Maxwell-Boltzmann and Fermi-Dirac statistics, in these cases the kernel n[z] 
falls off fast enough to make its integral finite. This is not the case for Bose-Einstein statistics for particles that 
decoupled in LTE while ultrarelativistic. In this case the asymptotic behavior of n[z] oc 1/z leads to an infrared 
enhancement as a consequence of the long-range of the kernel. In this case we keep the upper limit in eq. (j5.3p and 
carry out the integral in z, leading to [50| 



7" I 

"[^I'^^^CM.o 



dy 



ev -I 



1 — cos(^ uy) 



In [7 u e^ 
2C(3) 



(5.12) 



where C ~ 0.577216 ... is the Euler-Mascheroni constant as well as in the coefficient of the derivative term in eq. 
leading to 



/■7 " 

/ z n[z] dz 
Jo 



7 u — 'oo 7 ^ 



2C(3) 



(5.13) 



This infrared enhancement reflects the infrared divergence of the free streaming wavevector in Minkowski space-time 
found in ref. 37]. The argument of the logarithm in eq. (|5.12p clearly reveals that it is the cosmological expansion that 
yields an infrared cutoff. Taking u '--^ 1 [neglecting terms of 0(1/7)] we find a sliding wavevector at small scales for 
the Bose-Einstein case, namely 



K^(t,^) ~ 0.00424. . . v/W + C (l^j ■ j-^ [kpc]-i BE 



gd\3 m 



(5.14) 



The larger K^itf^q) leads to shorter free streaming lengths and to more power at small scales because free streaming 
smoothes out on shorter scales. Therefore, at small scales BE particles that decoupled while relativistic behave as 
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CDM. This conclusion is in agreement with the results in refs. |37l . |39| and is borne out in the numerical solution 
displayed in fig. (fTU]) . 

In all cases K^iteq) controls the decay of the gravitational fluctuations <i>(fc, u) for small scales and long times. The 
logarithmic behavior and the consequent increase of Ki^ yields an explanation of the enhancement at small scales over 
the FD and MB cases depicted in fig. (jTU]). 



Comparing statistics and initial conditions: 

The analysis presented above clearly indicates the differences arising from statistics and initial conditions. It proves 
convenient to compare the results of the numerical solution of Gilbert's equations for all cases considered parametrized 
by the dimensionless ratio 7. This comparison is depicted in fig. [TUlfor a wide range of 7. The differences from statistics 
and initial conditions are clearly displayed in this figure, and are in complete agreement with the analysis presented 
above. 
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FIG. 10: Comparison of exact numerical solutions of eq. (|2.47l ) and (|2.6ip for In |r(fc)| vs. 7 for Maxwell-Boltzmann particles, 
fermion and boson thermal relics with Gilbert and temperature initial conditions eas. p.65p and (|2.66p . 



For thermal relics that decoupled while relativistic, either FD or BE, the range of scales relevant for structure 
formation in which the linearized approximation is reliable corresponds to < 7 < 6, whereas for WIMPs, the 
relevant range of scales corresponds to 7 < 10~^. We have confirmed that in the range of cosmological relevance for 
all species, the second order approximation (|2.111[) for T{k) is very accurate and indistinguishable from the exact 
numerical solution in the range < 7 < 6. 

As shown in fig. [Tn]T(fc) features zeroes for MB and FD statistics for temperature initial conditions. However the 
value of 7 for FD at which T(k) vanishes corresponds to a sub-galactic scale, and for MB it corresponds to a sub-parsec 
scale, in either case well outside the regime of reliability of the linearized approximation, hence these features are not 
relevant for structure formation. Nevertheless, this figure highlights the main aspects discussed above: for a fixed 
value of 7 BE statistics favors the small momentum region, leads to a longer range memory kernel that falls off with 
a power law and yields the largest T{k) for a fixed value of 7, followed by FD statistics with a (slower) power law fall 
off and finally the MB distribution with an exponential fall off of the memory kernel and the smallest T{k) for fixed 
7- 



VI. CONCLUSIONS 



In this article we studied the evolution of gravitational and DM density perturbations from the coUisionless 
Boltzmann-Vlasov during matter domination, and obtained an exact expression for the transfer function T{k) for 
arbitrary distribution function of the decoupled DM particle and initial conditions. 

We have transformed the non-relativistic Boltzmann-Vlasov equation into an integro-differential equation which 
features a non-local kernel that describes the memory of gravitational clustering and yields corrections to the fiuid 
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description. This formulation lends itself to a systematic Fredholm expansion for the evolution of DM density and 
gravitational perturbations and T(k), and makes explicit the influence of the distribution function on T(k). 

Distribution functions that favor the small momentum region lead to longer range kernels, a persistence of the 
memory of the initial conditions and gravitational clustering resulting in an enhancement of T{k) at small scales. 

The natural scale of suppression of T{k) is determined by the free-streaming wave vector at matter-radiation 
equality, kfs{teq) = \/6/lfs{0) where the comoving free streaming wavevector is 



kfsit) = 



47rpoA/a(t) 



((I)' 



(6.1) 



the angular brackets refer to the average with the distribution function of the decoupled particle, and lfs{0) is the 
comoving free streaming distance traveled by the decoupled particle from the time of matter-radiation equality tgq 
until today. 



We find 
with 

0.325 



fc,,(te,)-^^ (6.2) 
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0.157 (^) ' (^) [kpc]-i FD thermal relic (^.3) 



0.175 ff-) ' (^) [kpc]-i BE thermal relic 



where gd is the number of relativistic species at decoupling. 

We provided a detailed numerical study for thermal relics, WIMPs and fermionic and bosonic particles that decou- 
pled while relativistically. The result of the numerical study of T{k) as a function of 7 = \/2k/kfs{teq) for the three 
cases and different initial conditions is displayed in fig. [TUl 

This study reveals that the first two terms in the Fredholm expansion yield a remarkable accurate approximation 
to T{k) in the range of scales of cosmological relevance for structure formation. In all cases considered we find that 
the exact solution differs from the second order approximation in less than 5% in the range < 7 < 6 corresponding 
to scales down to a fraction of the free-streaming length and much less than this on large scales for < 7 < 1. 

A simple and accurate approximation to T(k) is given by eqn. (|2.11ip . The second term in this expression includes 
corrections beyond the fluid approximation and higher order moments in the Boltzmann hierarchy and describes 
the memory of initial conditions and gravitational clustering. It explicitly depends on the distribution function of 
the decoupled particle. Non-local kernels with longer range lead to an enhancement of T{k) at small scales. FD 
and BE thermal relics feature kernels with longer range than WIMPs, with BE statistics (for relics that decoupled 
relativistically) leading to longer range and more power at small scales. 

This behavior is clearly exhibited in fig. [TOl 

For long scales k <^ kfs{teq) we find the behavior 




(6.4) 

where 70 is given by (|3.2ip for the initial conditions considered here independent of the statistics for thermal relics. 
We provide fits of T(fc), within a wide range of (small) scales for k > kfs{teq) in the intervals where it exhibits 
a simple power-like or exponential behavior for MB, FD and BE statistics: see eas. (|3.22p . (|4.13l) . (|4.14p . (|4.15p and 

Km . 

Although we do not attempt to provide a global fit with one function, it is clear that in the small scale region, 
the functional forms of T(fc) found in this article by exact numerical solution of the Boltzmann- Vlasov equation (in 
Gilbert's form) are very different from the fits by Bardeen eta/.J5l| often quoted in the literature. 

An important consequence of this study is a distinct imprint of the particle statistics and its distribution function 
on the transfer function at small scales, A <^ ^/s(0): a distribution function that enhances the small momentum 
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region yields a longer range memory of the initial conditions and gravitational clustering and an enhancement of the 
transfer function at small scales as depicted in fig. 1101 This result may prove important in the elucidation of the 
small scale structure of DM halos and perhaps lead to an explanation of the filamentary structures found in numerical 
simulations in ref. 52]. 

The tools provided in this article to study the transfer function for arbitrary distribution functions and general 
initial conditions, in particular the simple and remarkably accurate approximation to T{k) given by eqn. ()2.111|) 
allows a systematic and robust assessment of different DM candidates. 
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